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ABSTRACT 

An approximate Riemann solver for the equations of relativistic magnetohydrodynam- 
ics (RMHD) is derived. The HLLC solver, originally developed by Toro, Spruce and 
Spears, generalizes the algorithm described in a previous paper (Mignone & Bodo 
2004) to the case where magnetic fields are present. The solution to the Riemann 
problem is approximated by two constant states bounded by two fast shocks and sep- 
arated by a tangential wave. The scheme is Jacobian-free, in the sense that it avoids 
the expensive characteristic decomposition of the RMHD equations and it improves 
over the HLL scheme by restoring the missing contact wave. 

Multidimensional integration proceeds via the single step, corner transport upwind 
(CTU) method of Colella, combined with the contrained tranport (CT) algorithm to 
preserve divergence-free magnetic fields. The resulting numerical scheme is simple to 
implement, efficient and suitable for a general equation of state. The robustness of the 
new algorithm is validated against one and two dimensional numerical test problems. 
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1 INTRODUCTION 



Strong evidence nowadays supports the general idea that 
relativistic plasmas may be closely related with most of the 
violent phenomena observed in astrophysics. Most of these 
scenarios are commonly believed to involve strongly mag- 
netized plasmas around compact objects. Accretion onto 
super-massive black holes, for example, is invoked as the pri- 
mary mechanism to power hig hly energetic phenomena ob- 
served in active galactic nuclei, jMacch ettcll999l : lElvis et alJ 
l2002l:lMcKinnevl2005HShapircl2005l) . In this respect, the for- 
mation and propagation of relativistic jets and the accretion 
flow dynamics pose some of the most challenging and inter- 
esting quests in modern theoretical astrophysics. Likewise, a 
great deal of attention has been addressed, in the last years, 
to the darkling problem of gamma ray bursts (see for ex- 



ample |Meszaros & Rees! 


1994lMacFadyen & Wooslevlll999l 


iKonid & GranotJl2002l: 


Rosswoe et al.l|20od). whose mod- 



els often appeal to stro ngly relativistic coUimated outflows 
iAlov et al]l200ol . |2002|). Other attractive ex amples include 
pulsar wind nebulae jBucc iantini et al||2005h . microquasars 
i Meieil liool IMcKinnev fc Glmrlmiel 120041) . X-rav binaries 



Varniere et alJl2002t) an d stellar core collapse in the con- 



text of general relativity tBruennlll985t IDimmelmeier et alJ 
20021) . 

Theoretical investigations based on direct numerical 
* E-mail: mignone@to.astro.it 



simulations have paved a way towards a better understand- 
ing of the rich phenomenology of relativistic magnetized 
plasmas. Part of this accomplishment owes to the success- 
ful generalization of existing shock-capturing Godunov-type 
codes to relativistic magnetohydrodynamics (RM HD) (sec 
lKomissaroTlll999l iBalsardlioOll iDel Zanna et al.ll2003 . and 
reference therein). Implementation of such codes is based 
on a conservative formulation which requires an exact or 
approximate solution to the Riemann problem, i.e., t he de- 
cay o f a discontinuity separating two constant states iTord 
Il997l) . In terms of computational cost, employment of exact 
relativistic Riemann solvers may become prohibitive due to 
the high degree of intrinsic nonlinearity present in the equa- 
tions. This has focused most computational efforts towards 
the development of approximate solvers which, nevertheless, 
requi re knowledge of the e xact solution, at least on some 
level dMartf fc MiirleifeOOotl . The presence of magnetic fields 
further entangles the solution, since the number of decay- 
ing waves increas es from three to seven iAnile fc Pennisil 
Il987t lA~nile II1989T) . An exact analytical approach to the so- 
lution (which does not allow compound waves) ha s been 
recently prese nted in iGiacomazzo fc Rezzollal J2005I) . while 
iRomero et al.l i2005l) derived a special case where the veloc- 
ity and magnetic field are orthogonal. 

The trade-off between efficiency, accuracy and robust- 
ness of such approximate methods is still a matter of 
research. So lvers based o n loc al linearization ha ve been 
presented in iKomissarovl lll999l) (KO henceforth) , iBalsaral 
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fcOOlft (BA henceforth) and iKoldoba et"afl ©02). De- 
spite the higher accuracy in reproducing the full wave 
structure, these solvers rely on rather expensive charac- 
teristic decompositions of the Jacobian matrix. Conversely, 
the chara cteristic-free f o rmula tion of Harten-Lax-van Leer 
(HLL) of lHarten et alJ l)l983F) has gained increasing pop- 
ularity due to its ease of implementation and robustness. 
The HLL approach has been successfully applied to the 
RMHD equations bv iDel Zanna et aljEoM (dZBL hence- 
forth) as well as to the general relativistic c ase (see for ex- 
ample iGammie et aljl2003l : iDuez et"al . 2005) and to the in- 
vestigation of extragalactic iets. see lLeismann et al.l l)2005Fl . 

Besides the computational efficiency, however, the HLL 
formulation averages the full solution to the Riemann prob- 
lem into a single state, and thus lacks the ability to re- 
solve single intermediat e waves such as Alf v en, co ntact and 
slow discontinuities. In Mi gnone fc Bodo I (2005) (paper I 
henceforth) we proposed an approach that cured this defi- 
ciency by restoring the missing contact wave. The resulting 
sch eme genera li zed th e HLLC approximate Riemann solver 
by iToro et al.l il994T) to the equations of relativistic hy- 
drodynamics without magnetic fields. Here, along the same 
lines, we propose an extension of the HLLC solver to the rel- 
ativistic magnetized case. Similar w ork has been p rese nted 
in the context of classical MHD by iGurskil j2004) and iLi I 
(2005). 

The new HLLC Riemann solver is implemented in the 
framework of the corner transport upwind (CTU) method of 
IColellal jl99dl. coupled with the c onstrained transport (CT) 
evolution jEvans fc Hawlevlll988h of magnetic field. The al- 
gorithm naturally preserves the divergence-free condition to 
machine accuracy and is stable up to Courant number of 1. 

The paper is organized as follows. The relevant equa- 
tions are given in In J^we derive the new HLLC Riemann 
solver. Numerical tests, together with the implementation of 
the CTU-CT method are shown in 



2 THE RMHD EQUATIONS 

The motion of an ideal relativistic magnetized fluid is de- 
scribed by conservation of mass, 



d a (pu a ) = , 
energy-momentum , 



d 



(ph + \b\ z )u a v! 3 — b a b p + pr\ 



and by Maxwell's equations: 
d a {u a b -vPb a ) = . 



(1) 

= , (2) 
(3) 

see, for example. lAnile fc Pennisi I il987f) or lAnile I il988t) . 
In equations Q, and © we have introduced the rest 
mass density of the fluid p, the four velocity u a , the co- 
variant magnetic field b a and the relativistic specific en- 
thalpy h. The total pressure p results from the sum of ther- 
mal (gas) pressure p g and magnetic pressure |b| 2 /2, i.e., 
p = p g + \b\ 2 /2. In what follows we assume a flat metric, 
so that rf 10 — diag( — 1, 1, 1, 1) is the Minkowski metric ten- 
sor. Greek indexes run from to 3 and are customary for 
covariant expressions involving four- vectors. Latin indexes 
(from 1 to 3) describe three-dimensional vectors and are 
used indifferently as subscripts or superscripts. 



The four-vectors u a and b a are related to the spatial 
components of the velocity v = (v x ,v y ,v z ) and laboratory 
magnetic field B = (B x , B y , B z ) through 

u<* = 7 (l, v) , 

b a = jfvB, ^L+v(v-B] 

with the normalizations 
u a u a = —1 , u a b a = , 
IB' 2 



b a b a 



+ {v-B) 



(4) 

(5) 
(6) 



where 7 = (1 — v ■ v)^ 1 ^ 2 is the Lorentz factor. We follow 
the same conventions used in paper I, where velocities are 
given in units of the speed of light. 

Writing the spatial and temporal components of equa- 
tion @ in terms of the laboratory magnetic field yields 

'"' B = V x (v x B) , (7) 



Of 

V-B = 0, 



(8) 



i.e., they reduce to the familiar induction equation and the 
solenoidal condition. 

For computational purposes, equations are more 

conveniently put in the standard conservation form 



(9) 



dU ^ dF\U) _ 

dt 2^, Qx k 

k 

together with the divergence-free constraint ijHJ , where U = 
(D, mi , rn y , ra z , B x , B y , B z , E) is the vector of conservative 
variables and F k are the fluxes along the x k = (x, y, z) di- 
rections. The components of U are, respectively, the labo- 
ratory density D, the three components of momentum rrik 
and magnetic field Bk and the total energy density E. From 
equations Q, @ and the definitions Q one nas 



D 

m k 

E 
and 



= PI 1 

= {ph 1 2 +B 2 )v k -{v-B)B k , 



, 2 B 
= phy - Pg + — 



Dv x 



v 2 B 2 -(v-B) 2 



(10) 
(11) 
(12) 



F X {U) = 



rn x v x 



B x 



m y v x 



m z v x 



B,v x 



B u 
B x 



B x Vy 

B x v, 



V 



(13) 



Similar expressions hold for F V (U) and F Z (U) by cyclic 
permutations of the indexes. Notice that the fluxes entering 
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in the induction equation are the components of the elec- 
tric field which, in the infinite conductivity approximation, 
becomes 



n = -v x b 



(14) 



The non-magnetic case is recovered by letting B —* in the 
previous expressions. 

Finally, proper closure is provided by specifying an ad- 
ditional equation of state. Throughout the following we will 
assume a constant T-law, with specific enthalpy given by 



h = 1 + 



Pg 



(15) 



r-i P - 

where V is the constant specific heat ratio. 

2.1 Recovering primitive variables 

Godunov-type codes are based on a conservative formulation 
where laboratory density, momentum, energy and magnetic 
fields are evolved in time. On the other hand, primitive vari- 
ables, V — (p,v,p g , B), are required when computing the 
fluxes 11311 and more convenient for interpolation purposes. 

Recovering V from U is not a straightforward task in 
RMHD and different approaches have been suggested by pre- 
vious authors: BA used an iterative scheme based on a 5 x 5 
Jacobian sub- block of the system ©; KO solves a 3 x 3 non- 
linear s ystem of equations; dZ BL (the same approach is also 
used in lLeismann et all <|2005F0 further reduced the problem 
to a 2 x 2 system of nonlinear equations. Here we reduce this 
task to the solution of a single nonlinear equation, by prop- 
erly choosing the independent variable. If one sets, in fact, 
W — p/17 2 , S — m ■ B , the following two relations hold: 



E = W 



Pa 



1 - 



1 

5y2 



\B - 



S 2 
2W 2 



(16) 



\m\ 2 = (W+ \B\ 2 ) 2 (l - - |1 (2W + \B\ 2 ) . (17) 

Since at the beginning of each time step m, B and S 
are known quantities, equation 1171 allows one to express 
the Lorentz factor 7 as a function of W alone: 

( S 2 (2W+\B\ 2 ) + \m\ 2 W 2 y i 

Using the equation of state 1151 . the thermal pressure 
p g is also a function of W: 



P 9 {W) 



IY7 2 



(19) 



where F r = T/(T — 1) and 7 is now given by GHJ. Thus the 
only unknown appearing in equation (llfciP is W and 



f{W) = W -p 3 +[l-—) \B 



S 2 
2W 2 



E = (20) 



can be solved by any standard root finding algorithm. Al- 
though both the secant and Newton- Raphson methods have 
been implemented in our numerical code, we found the latter 
to be more robust and computationally efficient and it will 
be our method of choice. The expression for the derivative 
needed in the Newton scheme is computed as follows: 



df{W) _ dp a \B\ 2 d 1 S 2 
dW dW 7 3 dW W 3 



(21) 



where dp g /dW is computed from (1191 . whereas d'y/dW is 
computed from eq. 1181 : 



c?7 
dW 



dp g _ 7(1 + Ddj/dW) - 2Wdj/dW 

dw ~ r r 7 :i 



3 2S 2 {?>W 2 + 3W\B\ 2 + \B\ 4 ) + \m\ 2 W 3 
7 2W 3 (W+|B| 2 ) 3 



(22) 



Once W has been computed to some accuracy, the 
Lorentz factor can be easily found from 1181 . thermal pres- 
sure from 1191 and velocities are found by inverting equation 
CO}: 



Vk = 



W + 



1 ( s n 



(23) 



Finally, equation 11U1 is used to determine the proper 
density p. 



2.2 The Riemann Problem in RMHD 

In the standard Godunov-type formalism, numerical integra- 
tion of @ depends on the computation of numerical fluxes 
at zone interfaces. This task is accomplished by the (exact 
or approximate) solution of the initial value problem: 



U 



U(x,0) 



where U , 



if 
if 



x < x. 



X > X, 



(24) 



and U R i+ i are assumed to be piece- wise 
constant left and right states at zone interface i + |. The 
evolution of the discontinuity 1241 constitutes the Riemann 
problem. 

As in classical MHD, evolution in a given direction is 
governed by seven equations in seven independent conserved 
variables. Integration along the a;-direction, for example, 
leaves B x unchanged since the corresponding flux is identi- 
cally zero, eq. 1131 . The solution to the initial value problem 
1241 results, therefore, in the formation of seven waves: two 
pairs of magneto-acoustic waves, two Alfven waves and a 
contact discontinuity. 

The complete analytical solution to the relativistic 
MHD R iemann problem has been recen tly derived in closed 
form bv lGiacomazzo fc Rezzolial (I2005T) . A number of prop- 
erties regarding s i mple waves are also well e stablished, see 
I Anile fc Pennisi I (Il987f) and lAnile I dl989l) . iRomero et all 
l2005ir "discuss the case in which the magnetic field of the 
initial states is tangential to the discontinuity and orthogo- 
nal to the flow velocity. 

General guidelines, relevant to the present work, follow 
below. Across a magneto-acoustic (fast or slow) shock, all 
components of V can change discontinuously. Thermody- 
namic quantities (e.g., p and p g ) are continuous through a 
relativistic Alfven wave (as in the classical case), but con- 
trary to the classical counterpart, the magnetic field is ellip- 
tically polarized and the normal com ponent of the velocity 
is discontinuous llKomissarovl I199"t1) . Through the contact 
mode, only density exhibits a jump while thermal pressure, 
velocity and magnetic field are continuous. 
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For the special case in which the component of the mag- 
netic field normal to a zone interface vanishes, a degeneracy 
occurs where tangential, Alfven and slow waves all propa- 
gate at the speed of the fluid and the solution simplifies to a 
three-wave pattern. Under this condition, the approximate 
solution outlined in paper I can still be a pplied with mi- 
nor m odifications, see H3.2l in this paper and lMignone et all 
(2005). 



F* L X R {X* - Xl) + F R \ L {X R - A*) 



3 THE HLLC SOLVER 

The derivation of the HLL and HLLC approximate Riemann 
solvers has already been discussed in paper 1 and will not 
be repeated hereafter. 

Following the same notations, we approximate the so- 
lution to the initial value problem 1241 with two constant 
states, U* L and U* R , bounded by two fast shocks and a con- 
tact discontinuity in the middle. We write the solution on 
the x/t = axis as 



1/(0, t) = < 



U* L 
Ur 

Ur 



if 

if 

if 
if 



Al ^ , 
A L < 5? A* 
A* < s= Xr 

Xr^O, 



(25) 



where Xl and Xr are, respectively, the minimum and max- 
imum characteristic signal velocities and A* is the velocity 
of the middle contact wave. The corresponding inter-cell nu- 
merical fluxes are: 



f=< 



F L 
Ft 
Fr 
Fr 



if 
if 
if 

if 



X L > , 
X L < < A* 
A* < < X R 

Xr < . 



(26) 



The intermediate fluxes F* L and F* R are expressed in 
terms of U* L and U* R through the Rankine-Hugoniot jump 
conditions: 



Xl (Ul 


-Ul) 


= Ft 


-Fl , 




A* (U* R 


-ut) 


= F* R 


-Ft , 


(27) 


X R {Ur 


-u R ) 


= F R 


-Ft , 





where, in general, F* L>R + F{U* L , R )- 

The consistency condition is obtained by adding the 
previous equations together: 



(A* — Xl)U*l + {Xr — X*)U* R = jjhii 
Xr — Xl 

where 

hu XrUr - XlUl + F L - F R 



U 



Xr — Xi 



(28) 



(29) 



is the state integral average of the solution to the Riemann 
problem. 

Similarly, if one divides each expression in eq. (1271 by 
the corresponding A's on the left hand sides and adds the 
resulting expressions, 



Xi 



Xr 



= X* F 



with 



phii _ X R F l — XlFr + X r Xl{Ur — Ul) 
Xr — Xl 



(30) 



(31) 



being the flux integral average of the solution to the Rie- 
mann problem. 

Since the sets of jump conditions across the contact 
discontinuity differ depending on whether B x vanishes or 
not, we proceed by separately discussing the two cases. In 
either case, the speed of the contact wave is assumed to be 
equal to the (average) normal velocity over the Riemann 
fan, i.e. A* = v%. The normal component of magnetic field, 
B x , is assumed to be continuous at the interface, so that 
Bx = B X} l = B x< r can be regarded as a parameter in the 
solution. 



3.1 Case B * 

We start by noticing that equations 12811 and 1301 1 provide 
a total of 14 relations. Six additional conditions come by 
imposing continuity of total pressure, velocity and magnetic 
field components across the contact discontinuity. This gives 
us a freedom of 20 independent unknowns, 10 per state; we 
choose to introduce the following set of unknowns for each 
state 



, Vy, By, B* 



i*, E*, p*} 



(32) 



The normal component of momentum {m%) is not an inde- 
pendent variable since we assume, for consistency, that 



m* = {E* +p*)v* - {v* ■ B*)B* 



(33) 



The previous relation obviously holds between conservative 
and primitive physical quantities. We point out that the 
choice 13211 is not unique and alternative sets of independent 
variables may be adopted. 

According to the previous definitions, the state vector 
solution to the Riemann problem is written as 



,m y ,m* z ,B*,Bl,E* 



U* = \D*,ml, 

while the flux vector, eq. 1131 , becomes 
/ D*v* 
. . BtBl 



(34) 



F* 



(7* 



B*v* {v* ■ B*) +p* 



B%v* {v ■ B ) 



(T) 2 
BIB* Z 
(7*) 2 



B*v* 



B*v*z («* 



B'Vy 



V 



/ 



(35) 



As in paper I, we adopt the convention that quantities with- 
out the L or R suffix refer indifferently to the left (L) or right 
{R) state. 

The six conditions across the contact discontinuity are 
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y y,L 



(36) 

Pl = Pr ■ 

For these quantities the suffix L or R is thus unnecessary. 

From the transverse components of the magnetic field in 
the state consistency condition (J2HJ, one immediately finds 
that 



B,, 



Bt — B, 



(37) 



Thus the transverse components the magnetic field are given 
by the HLL single state. Similarly, from the fifth and sixth 
components of the flux consistency condition 1301 one can 
express the transverse velocity through 



B X Vy By V x 



g 



B*v* = B*v* 



(38) 



where Fg l * and Fg 1 ' are the B y - and B z - components of the 
HLL flux, eq. 1311 . Simple manipulations of the normal mo- 
mentum and energy components in equation 1281 together 
with 1331 yield the following simple expression: 



E v x +p v x 



B x 



(v* -B*) => 



(39) 



Similar algebra on the momentum and energy compo- 
nents of the flux consistency condition 1301 leads to 



Bl (v*-B*] 



1 J* 

7* 



■P*-F%1 =0 



(40) 



where l/( 7 *) 2 = 1 - (v*) 2 - (v* y f - «) 2 . 

Now, if one multiplies equation 14011 by v* and subtracts 
equation l|39[l . the following quadratic equation may be ob- 
tained: 



a(v x ) 2 + bv* + c ■■ 
with coefficients 
a = F|" - 1 

b = -F™1- 
—Ml 



. 



phll 
B ± > 

B] 

hll 



(41) 



(42) 



In the previous expressions B Ml = (0,B Ml ,B Ml ), F Ml ± = 
(0, Fg , Fg" ) . Similar arguments to those presented in paper 
I lead to the conclusion that only the root with the minus 
sign is physically admissible. 

Once v x is known, Vy and v* are readily obtained from 
13811 , p* is computed from 1401 , while density, transverse mo- 
menta and energy are obtained using the Rankine-Hugoniot 
jump conditions across each fast wave: 

A 



D* = 



A - 



-D 



-BZ 



+ {v* 



B*)v* y 



Xm y — F n 



-BZ 



BZ 

W) 



7 + (v* ■ B*)v* 



Am, 



F„. 



E* 



\-v% 

\E — m x + p*v* — (v* ■ B*)B X 

x — vt 



(43) 



(44) 



(45) 



(46) 



In equations 1441 and 1451 . F my and F mz are, respectively, 
the m y - and m z - components of the flux, eq. 1131 . evaluated 
at the left or right state. As in paper I, we have omitted the 
suffix L or R for clarity of exposition. 



3.2 Case B* = 

For vanishing normal component of the magnetic field a de- 
generacy occurs where the Alfven waves and the two slow 
magnetosonic waves propagate at the speed of the contact 
discontinuity. For this case the approximate character of the 
HLLC solver offers a better representation of the exact so- 
lution, since the Riemann fan is comprised of three waves 
only. At the contact discontinuity, however, only the nor- 
mal component of the velocity v x and the total pressure 
p are continuous (KO). The remaining variables experience 
jumps. This only adds 2 constraints to the 14 jump condi- 
tions, leaving a freedom of 8 unknowns per state. However, 
the transverse velocities v y and v z do not enter explicitly 
in the fluxes 1351 and the jump conditions can be written 
entirely in terms of {D*,v x ,m y ,rn%,B*,B*,E*,p*}, i.e. 8 
unknowns per state. Straightforward algebra shows that the 
coefficients of the quadratic equation 1411 are now given by 



a = F F 



b = -F n 



E n 



(47) 



i.e., they coincide with the expressions derived in paper I. 
The root with the minus sign still represents the correct 
physical solution. Once v x is found, the total pressure p* is 
derived from 



? hll 



p = -F E v x + F m 
and the normal momentum 1331 becomes 
m* = (E* +p*)v* . 



(48) 



(49) 



The remaining quantities are easily obtained from the 
jump conditions: 

A — u~ 



D 



-D 



E" 



X-vi 



X-v x 

XE - m x +p*v* 

X — Vx 
X — v x 



(50) 
(51) 
(52) 
(53) 



3.3 Remarks 

The expressions derived separately in H3.1l and H3.2l are suit- 
able in the B x ^ and B x —* cases, respectively. Although 
other degeneracies may be present (see KO for a thorough 
discussion) no other modifications are necessary to the algo- 
rithm. Before testing the new solver, however, a few remarks 
are worth of notice: 

(1) The solutions derived separately for B x 7^ and the 
special case B x = automatically satisfy the consistency 
conditions 1281 and 1301 by construction; 

(2) In the limit of zero magnetic field, the expressions 
derived in H3.2l reduce to those found in paper I; 

(3) In the classical limit, our derivation does not coin- 
cide wi t h the app roximate R iemann solvers constructed by 
iGurskil <2004t) or ILTI l|2005l) . The reas on for this dis crep- 
ancy stem s from the fact that both IGurskil (120041) and 

ED 

(2005) assume that transverse momenta and velocities 
are tied by the relation m y>z = p*Vy tZ . Although certainly 
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true in the exact solution, this assumption reduces, in the 
HLLC approximate formalism, the number of unknowns 
from 10 to 8 (when B x =^ 0) thus leaving the systems of 
jump conditions 1271 overdetermined. Should this be the 
case, the number of equations exceeds the number of un- 
knowns and the integral relations across the Riemann fan in- 
evitably break down. This explains the inconsistencies found 
in Li's and Gurski ' s der ivations and further discussed in 
iMivoshi fc Kusanol J2005t) . 

Therefore, in the classical limit, our expressions automati- 
cally imply rriy iZ ^ p*Vy iZ and the correct expressions for the 
transverse velocities are still given by 11381 . whereas trans- 
verse momenta should be derived from the jump conditions 
accordingly. Furthermore, contrary to Li's misconception, 
consistency with the jump conditions requires that the mag- 
netic field components be uniquely determined by 13711 and 
no other choices are thus possible. 

(4) The reader might have noticed that in the limit of 
vanishing B x , some of the expressions given in H3.1l do not 
reduce to the those found in H3.2 | This prop erty also pe rsists 
in the classical limit, see lGurs ki (2004). and^T] {2005). The 
reason for this discrepancy relies on the assumption of conti- 
nuity of the transverse components of magnetic field across 
the tangential wave A*: when B x — > 0, a degeneracy occurs 
where the tangential, Alfven and slow waves all propagate at 
the speed of the fluid and the solution simplifies to a three- 
wave pattern. In the exact solution, the continuity of B y and 
B z across the tangential wave is lost since the middle state 
bounded by the two slow waves becomes singular. 

(5) Lastly, we note that in both the classical and relativis- 
tic case the transverse velocities given by eq. 1381 become 
ill-defined as B x — » 0. However, in the classical case, the 
terms involving Vy or v z in the flux definitions remain finite 
as B x — > 0. Conversely, this is not the case in RMHD for 
arbitrary orientation of the magnetic field as one can see, 
for example, using eq. 1441 : 



(B h z 



jphll\/ jphll qHII 

B z ) V By z 



TT'hll T)hll\ 



B X (X - «*) 



+ 0(1) (54) 



as B x — > 0. Fortunately, for strictly two dimensional flows 
(e.g. when B z — v z = 0) the leading order term vanishes 
and the singularity is avoided. In the general case, however, 
we conclude that more sophisticated solvers should allow 
the presence of rotational discontinuities in the solution to 
the Riemann probl e m. Th is has been done, for example, by 
IMivoshi fc Kusanol (120051) in the context of classical MHD. 



3.4 Wave Speed Estimate 

The full characteristic decomposition of the RMHD equation 
(i.e. the eigenvalues and eigenvectors of t he Jacobian ma- 
trix d F x /d U) was extens ively analyzed bv lAnile fc Pennisi I 
( 1987) and lAnile I l)l989h . In the one-dimensional case the 
Jacobian matrix can be decomposed into seven eigenvectors 
associated with four magnetosonic waves (fast and slow dis- 
turbances), two Alfven waves and one entropy wave propa- 
gating at the fluid velocity. The eigenstructure is therefore 
similar to the classical case and it can be shown that the or- 
dering of the v arious speed s and corresponding degeneracies 
are preserved jAnile Ill989h . 

Since the HLLC approximate Riemann solver requires 



an estimate of the outermost waves, the right and left-going 
fast shock speeds identif y the n ecessary characteristic veloc- 
ities. Thus we set (IDavislllfl 



(55) 



X L =mm(\-(V L ),\-{V R )) , 

\ R = max(\ + {VL),\+(V R )) , 

where A_ and A+ are the minimum and maximum roots of 
the quartic equation 

ph(l - c 2 s )a 4 = (1 - A 2 ) [(|6| 2 + phc 2 s )a 2 ~ c 2 B 2 ] , (56) 

with a = 7(A — v x ), B = b x — \b°. In absence of mag- 
netic field, both the (left and right-going) slow and fast 
shocks propagate at the same speed and equation 1561 re- 
duces to the quadratic equation (22) shown in paper I. When 
B 7^ 0, no simple analytical expression is available and solv- 
ing 1561 requires nume rical or rather cumbers ome analytical 
approaches. Recently, iLeismann et alJ ll2005l) proposed ap- 
proximate simple lower and upper bounds to the required 
eigenvalues. Here we choose to solve eq. (1561 by means of 
analytical methods, where the quartic is reduced to a cubic 
equation which is in turn solved by standard methods. 

There are special cases where it is possible to handle 
some of the degeneracies more efficiently using simple ana- 
lytical formulae: 

• for vanishing total velocity, equation 1561 reduces to a 
bi-quadratic, 



(ph+\b\ 2 )\ 4 



<\ 2 + phc 2 + B 2 c 2 )A 2 + c s Bl 







(57) 



• for vanishing normal component of the magnetic field, 
equation 1561 yields a quadratic equation: 



02A + ai A + ao = 



(58) 



with a 2 = ph[c 2 + 7 2 (l-c 2 )] + Q, ai = ~2ph-y 2 v x (l - c 2 ), 
a = pfe[-c 2 + 7 2 t; 2 (l-c 2 )] - Q and Q = \b\ 2 -c 2 (v ± -B ± ) 2 . 

For all other cases we solve the quartic equation 1561 . 



3.5 Positivity of the HLLC scheme 

The set of physically admissible conservative states, G, iden- 
tify all the LT's yielding positive thermal pressure p g and to- 
tal velocity \v\ < 1, according to the procedure outlined in 
>I2.1I Thus the positivity of the HLLC approximate Riemann 
solver requires that 

• both left and right intermediate states U* L and U* R be- 
long to G; 

• the first-order scheme yields updated conservative 
states that are in G. 

Unfortunately, the mathematical proof of the positivity 
of the HLLC scheme presents remarkable algebraic difficul- 
ties. In absence of the singular behavior described in < j3.3l 
investigations have been carried at the numerical level by 
verifying that each intermediate state U* correspond to a 
primitive, physically admissible state. In all the tests pre- 
sented in this paper and several others not discussed here, 
the scheme did not manifest any loss of positivity. However, 
in the general three-dimensional case when B x , B y ,B z ^ 0, 
the terms involving B x in the expressions for the transverse 
momenta may become arbitrarily large as B x — > and a loss 
of positivity can be experienced. 
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4 ALGORITHM VALIDATION 

4.1 Corner Transport Upwind for relativistic 
MHD 

The RMHD equations © are evolved in a conservative, di- 
mensionally unsplit fashion: 



U 



n+l 



(59) 



where the Cs are Godunov operators 



At / *,n+i 



f 



C?) 



(60) 
(61) 



At 

' Ay 3 \'i>3+i 

and U n is the set of volume-averaged conservative variables 

U n = (^D,m,B,E^j at time t = t n . Here B denotes 

the zone-averaged magnetic field. For clarity of exposition 
we will omit, throughout the following, integer-valued sub- 
scripts and retain only the half- integer notation to de- 
note zone edge values. 

The fluxes appearing in equations 1601 and I6H are 
computed by solving, at each zone interface, a Riemann 
problem with suitable time-centered left and right input 



states. For example, we obtain / 



i+i 



as the HLLC flux 



and V ? .respectively. 



with input states given by V . 

Computation of time-centered left and right zone edge 
va lues pro c eeds using the corner transport upwind (CTU) 
of IColellal dl990h. recently extende d to relativistic hydro- 
dv namics bv [ Mi gnone et alJ (l2005aT) and to classical MHD 
by iGardiner fc Stonel <l2005l) . Here we generalize the CTU 
approach to relativistic MHD by following a slightly differ- 
ent approach , although equivalent to the guidelines given in 
IColellal (Il99(j> . For the sake of conciseness, only the essential 
steps will be described hereaft e r. The unfamiliar r eader is re- 
ferred to the work of IColellal il99(J> . ISaltzmarJ Jl994l) and 
IGardiner fc Stone! i2005T) for more comprehensive deriva- 
tions. 

In our formulation, second-order accurate left and right 
states are sought in the form 



V 



~±- 



8 V V n 



2 ' 3±3> s 2 ' 

where we take S = L (S — R) with the plus (minus) sign. 
The slopes S X V" and S y V n are computed at the beginning 
of the time step using, for example, the monotonized central- 
difference (MC) limiter: 



5 x q n = s i min I 2|A<#], 
where q 6 V and 



2\Aq r l 



m+i - <ii 



Aq± = ± (q i± i - Qi ) , Si = 



sign(A(4)+sign(AgI! 



An alternati ve smoother pr escription is given by the har- 
monic mean (Ivan Leerlll977t) : 

2max(0, Aq+Aq-) 



similar manner. Additional forms of limiting may be adopted 
if necessary, see jATland gXj 

The cell- and time- centered values on the right hand 
sides of equations I62H are computed from a Taylor expan- 
sion of the conservative variables, i.e. 



U J 



U 



U n + —— = U n - — 

2 dt 2 



n At dU _ „ At 
2 dt 2 



OF dF y 



dx 



8F X 



+ 



+ 



OF" 



,(66) 



.(67) 



dx dy 

Following IColellal JT990V we approximate the spatial deriva- 
tive in the direction normal to a zone interface (denoted with 
a hat) with the Hancock step already introduced in paper I, 



d F\ F ' x { v : + ^)- FX { v ^ 

dx Axi 



(68) 



whereas the derivative in the tangential direction is com- 
puted in an upwind fashion using a Godunov operator: 

8F y „,,„ At 



At 



dy 



-C v ' n = 



The state U y ' n+ z is obtained by similar arguments by inter- 
changing the role of normal and tangential derivatives. We 
would like to point out that the Godunov operators used in 
the predictor step involve left and right states computed at 
t = t n (and not at t = t n+ i as in lGardiner fc Stonel I^PORTl ): 

v; ±i , s = v-±^-. (70) 



At ( y, n 

Ay, \ J i+h 



(69) 



This choice still makes the scheme second-order accurate in 
space and time and was found, in our experience, to yield 
a more robust algorithm. Besides, our CTU implementation 
does not require a primitive variable formulation, thus of- 
fering ease of implementation in the context of relativistic 
hydro and MHD, where the Jacobian dF/dU is particularly 
expensive to evaluate. 

Note that a total of four Riemann problems are involved 
in the single time step update 159H . It can be easily verified 
that for one-dimensional flows, the corner transport upwind 
method outlined above reduces to the scheme presented in 

(62) P a P er L 

Finally, the choice of the time step At i s based on the 
Courant-Friederichs-Lewy (CFL) condition dCourant et alJ 
1928): 



At = CFL x min 



Ax 



Ay 



i,j Vmax(|A||,|A-|)'max(|A||,|A^|) 



(71) 



(63) 



(64) 



(65) 



Aq+ + Aq- 

Equation l|63|l provides smaller dissipation at disconti- 
nuities, whereas equation 165H was found to give less oscil- 
latory results. Interpolation in the y-direction is done in a 



where < CFL < 1 is the Courant number and |Af, R \, 
r\ are ^ ne zone interface wave speeds computed in the x 
and y directions according to I55> . 



4-1.1 Contrained Transport Evolution of the Magnetic 
Field 

It is well known that multidimensional numerical schemes 
do not generally preserve the solenoidal condition, eq. JSJ, 
unless special discretization techniques are employed. In 
this respect, several approaches have been suggested in 
the context of the classical MHD equations (iTothl 119971 : 
lLondrillo fc Del Zann al l2000l) and some of them have been 
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recently extended to the relativistic case, see dZBL. Here 
we ad opt the constrained transpo rt (CT) ( lEvans fc Hawlevl 
T988!) and follow the approach of lBaIsara^^SplcerTi^9Sr i 
for its integration in Godunov-type schemes. 

In the CT approach a new staggered magnetic field vari- 
able is introduced. In this representation, the components of 
the magnetic field are treated as area-weighted averages on 
the zone faces to which they are orthogonal. Thus, B x is 
collocated at (i + whereas B y at + ~). No jump is 
allowed in the normal component of B at a zone boundary, 
consistently with the well posedness of the Riemann prob- 
lem presented in >12.2I and [0 Transverse components may 
be discontinuous. 

In this formulation, a discrete version of Stoke's theorem 
is used integrate the induction equation J7J. For example, 
after the predictor steps 166H and 167H . we update the face- 
centered magnetic field according to 



B n+ \ 



B" 



= B n -i ] 



= B n ., i 



+ 



At n 
2A~y~~ 

At n 
2Axi 



n z 



n z 



n z 



(72) 



and similarly after the corrector step. The electromotive 
force fi is collocated at cell corners and is computed by 
straightforward arithmetic averaging: 

Q i+h,3+h = Z ' (73) 



where, fl z 



- f x ' n , and Q z . , , 



f v '™ , are the 



z components of the electric fields available at grid interfaces 
during the upwind step. Despite its simplicity, eq. (1731 lacks 
of directional bias and more sophisticated algorithms may 
be used to incorporate upwind informa t ion in a consistent 
way, see lLondrillo fc Del Zannal J2004t) . iGardiner fc Stond 
(2005). For ease of implementation we will not discuss them 
here. 

It is a straightforward exercise to verify that the V ■ B = 
condition is preserved from one time step to the next one, 
due to perfect cancellation of terms. Notice also that, since 
B x is continuous at the (i + interface, only B y and B z 
need to be interpolated during the reconstruction procedure 
in the a;-direction. A similar argument applies to B x and B z 
when interpolating along the y coordinate. 

Since equation I59H evolves volume- averaged quantities, 
the zone-averaged magnetic field, B, is computed at the 
beginning of the time step from the face-averaged magnetic 
fields using linear interpolation: 



B x — 



By = 



B v,i+i + B y,j- 



(74) 



(75) 



Equations 1731 . 1741 and 1751 1 are second-order accurate in 
space. 

4-1-2 Summary 

We summarize our CTU constrained transport algorithm by 
the following steps: 

(1) At the beginning of the time step, form the volume 
averages 1741 and 1751 from the face centered magnetic field. 



(2) Compute x and y limited slopes by interpolating cell 
centered primitive variables according to eq. 1631 or 1651 . 

(3) Make a sweep along the x direction. Form left and 
right states using the first of eq. 17011 with B^ i+ i L = 
B™ i R equal to the x component of the face centered mag- 
netic field; 

- use the Hancock step l|68|l to compute the x derivative 
in eq. 1661 and add the resulting contribution to U x ' n+ ? ; 

- compute the C x ' n Godunov operator by solving Rie- 
mann problems at the (i + |, j) interfaces and add the 

resulting contribution to 

(4) Make a sweep along the y direction. Form left and 
right states using the second in eq. 170H with L = 
By j+i r equal to the y component of the face centered mag- 
netic field; 

- obtain the C y,n Godunov operator 1691 by solving 
Riemann problems at the + |) interfaces; add the 

resulting contribution to U*'™ + 5 . 

- use the Hancock step relative to the y direction to 

compute the y derivative and add it to !/| ''J 1 3 ; 

(5) Compute the time-centered area weighted magnetic 
field using Stoke's theorem 1721 . This concludes the predic- 
tor step. 

(6) Make a sweep along the x direction with left and 
right time-centered states given by the first equation in 1621 

* , = B™ +5 , „ equal to the time centered face- 



with B 

x 

averaged magnetic field computed via Stoke's theorem. Ob- 
tain the C x,n+ ^ Godunov operator. 

(7) Repeat the previous step by sweeping along the y di- 
rection. Compute the £ v,n+ 2 Godunov operator. 

(8) Update the cell-centered conservative variables using 
eq. 159H and the face-averaged magnetic field using Stoke's 
theorem. 



4.2 One-dimensional test problems 

One-dimensional problems are specifically designed to verify 
the ability of the algorithm in reproducing the exact wave 
pattern. In what follows we present four shock-tube tests, 
already introduced by BA and dZBL, with left and right 
states given in Table Q Computations are performed on the 
interval [0, 1] and the initial discontinuity is placed at x = 
0.5. The final integration time is t — 0.4. Note that the 
constrained transport algorithm is unnecessary, since eq. JSJ 
is trivially satisfied in one-dimensional flows. 



4-2.1 Problem 1 

The fi rst test problem, initially prop osed bv [ van Putter] 
(1993), is a relativistic extension of the iBrio fc Wu I ( ^8^1 
magnetic shock tube. In analogy with the classical case we 
use the ideal equation of state I15H with specific heat ratio 
r = 2. The breakup of the initial discontinuity sets up a left- 
going fast rarefaction wave, a left- going compound wave, a 
contact discontinuity, a right-going slow shock and a right- 
going fast rarefaction wave. 
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7 
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1 
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10 


-7 


-7 



Table 1. Initial conditions for the one-dimensional shock tube 
problems presented in the text. In all test problems we adopt 
a resolution of 1600 uniform computational zone, covering the 
interval [0, 1]. Integration is carried until t = 0.4. 



We compare, in Fig. Q the results obtained with the 
first-order HLL and HLLC solvers on 100 uniform com- 
putational zones. The exact solution (given by the solid 
line) was obtained using the n umerical code available from 
iGiacomazzo fc Rezzollal <l2005f) . The left going compound 
wave located at x ~ 0.5 is only visible in the numerical 
integration since the code used to generate the analytical so- 
lution (shown as the solid line in Fig.0 does not allow com- 
pound structures by construction. As expected, the HLLC 
Riemann solver attains sharper representation of the contact 
discontinuity when compared to the HLL scheme. Because 
of the reduced smearing in proximity of the contact wave, 
neighboring structures such as the compound wave on the 
left and the slow shock on the right can be better resolved 
when using the HLLC solver. Computations at different res- 
olutions show, in fact, that the L-l norm errors in density 
are reduced by roughly 20 -f- 30% (see left panel in Fig. 0, 
with Li(%) being, respectively, 0.53 and 0.74 for the HLLC 
and HLL solver at the highest resolution employed (6400 
zones). 

Fig.[3]shows the results obtained with the second-order 
scheme with the MC limiter, eq. (16311 . and the same Courant 
number, CFL = 0.8 on 1600 grid points. A direct compar- 
ison with the exact solution shows that all discontinuities 
are correctly captured and resolved on few computational 
zones, owing also to the presence of a compressive limiter. In 
this respect, our second-order HLLC scheme provides sim- 
ilar results to those obtained with the third-order central 
ENO-HLL scheme by dZBL. 

The L-l norm errors computed at different resolutions 
with the two different solvers differ by ~ 10 -f- 20%, see left 
panel in Fig. |1] When compared to the more sophisticated, 
characteristic-based algorithm presented in BA, our results 
show slightly sharper representation of the right-going slow 
shock and the contact discontinuity. Small overshoots ap- 
pear in the Lorentz factor profile at the left going compound 
wave and the right going slow shock. More diffusive slope 
limiters do not exhibit this feature. 

4.2.2 Problem 2 

The resulting wave pattern for this configuration is com- 
prised of two left-going rarefaction fans (fast and slow) and 
two right-going slow and fast shocks. The specific heat ratio 




0.0 I , , , i , , , i , , , i , , , i , , 

0.0 0.2 0.4 0.6 0.8 1.0 

x 

Figure 1. Comparison between the first-order HLL (dotted line) 
and the HLLC (dashed line) method for the first shock tube 
problem at t = 0.4. Only density profiles are shown. Compu- 
tations were performed on 100 computational zones with CFL 
= 0.8. The solid line gives the analytic solution as computed by 
IGiacomazzo fc R czzolla ( 20Q3). The major difference between the 
two approaches is the resolution of the contact wave. 





L,(%), PI 








10 


• HLL 






* HLLC 




1 



100 1000 100 1000 

Ax -1 Ax~' 

Figure 2. Discrete Ll-norm density errors (in percent) computed 
for the first-order scheme at different grid resolutions using the 
HLLC (asterisks) and HLL (filled circles) solvers. Computation 
have been performed for the first (left panel, PI) and second 
(right panel, P2) problems on 50, 100, 200, 400, 800, 1600, 3200 
and 6400 zones with CFL = 0.8. 

used for this calculation is F = 5/3. The weak slow rarefac- 
tion located at x ~ 0.53 and the slow shock at x ~ 0.86 are 
separated by a contact discontinuity where the proper den- 
sity changes by a factor of ~ 7. The velocity on either side 
of the contact wave is mildly relativistic, with a maximum 
Lorentz factor of ~ 1.36. 

The improvement offered by the HLLC Riemann solver 
over the HLL approach in the resolution of the contact wave 
is evident from Fig.|S] where we compare the density profiles 
obtained with the first order schemes against the analytical 
solution. 

Computations obtained with the second-order limiter 
<63> show excellent agreements with the analytical profiles, 
see Fig. |S| Our single-step HLLC scheme attain consider- 
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Figure 3. Relativisitc Brio-Wu shock tube problem. The second-order scheme with the HLLC Riemann solver on 1600 grid points and 
the MC limiter was used. From left to right and top to bottom: proper density, thermal pressure, Lorentz factor, normal and transverse 
velocity components and transverse magnetic field. The Courant number is 0.8. 
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Figure 4. Discrete Ll-norm error (10 2 ) for density computed 
for the second order scheme at different resolutions, see Fig. [5] 




ably sharper resolution than the results obtained by previ- 
ous calculations. The two right-going shocks, for instance, 
are smeared over ~ 3 grid points, approximately half of the 
resolution shown in BA and dZBL. Moreover, the smearing 
of the contact wave is considerably reduced when compared 
to the HLL scheme in dZBL (~ 10 zones vs. ~ 14). Similar 
overshoots, though, appear at the right of contact mode. 

The discrete L-l errors for different grid sizes are shown 
in the right panel of Fig. 3] where, at the maximum reso- 
lution employed (6400 zones) the HLLC and HLL errors 
reduce to 0.17% and 0.25%, respectively. 



Figure 5. Comparison between the first-order HLL (dotted line) 
and the HLLC (dashed line) method for the second shock tube at 
t = 0.4. Density profiles are shown. Computations were performed 
on 100 computational zones with CFL = 0.8. The solid line gives 
the a nalytic solution as computed by [Ciacomazzo fc RezzoIIal 
(2005). 



4.2.3 Problem 3 

The configuration for this test is similar to the previous 
problem, but a higher pressure jump separates the initial left 
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Figure 6. Solution of the mildly relativistic blast wave problem (test 2) computed with the second-order HLLC scheme and the MC 
limiter. A Courant number of 0.8 and 1600 grid zones were used in the computation. 



and right states, see Tabled Only the second-order scheme 
with the Van Leer limiter (1651 and a Courant number of 0.8 
has been employed. The ideal equation of state l|15p with 
r = 5/3 is used. The ensuing wave pattern shows a stronger 
relativistic configuration, with a maximum Lorentz factor 
of ~ 3.37, see Fig. [7j The presence of magnetic fields makes 
the problem even more challenging than its hydrodynam- 
ical counterpart (see test 3 in paper I), since the contact 
wave, slow and fast shocks now propagate extremely close 
to each other. As a result, a thin density shell sets up be- 
tween the contact mode and the slow shock. The higher 
compression factor (more than 100) follows from a more 
pronounced relativistic length contraction effect. At the res- 
olution of 1600 grid zones, the relative error in the density 
peak (pmax ~ 9.98) is 1.2%. A second thin shell-like struc- 
ture forms between the slow and fast shocks, as can be seen 
in the profiles in Fig. |7| The peaks achieved in the trans- 
verse components of velocity (« —0.37) and magnetic field 
(fa 8.95) achieve, respectively, 87% and 95% of their exact 
values. The small shell thickness, however, still prevents a 
clear resolution of the two right going shocks, visible in the 
exact solution. This demonstrates that relativistic magne- 
tized flows can develop rich and complex features difficult 
to resolve on a grid of fixed size. Similar conclusions have 
been drawn by previous investigators. 

Results obtained with the HLL solver (not shown here) 
indicates that the resolution attained at the contact discon- 
tinuity is equivalent. Therefore, as it was also pointed out 
in paper I, we conclude that, for strong blast waves where 
relativistic contraction effects produce closely moving dis- 
continuities, the HLL and HLLC schemes produce nearly 
identical results. 



4.2.4 Problem 4 

The collision of two relativistic streams is considered in the 
fourth test problem. The initial impact produces two strong 
relativistic fast shocks propagating symmetrically in oppo- 
site direction about the impact point, x = 0.5, see Fig. |S] 
Two slow shocks delimiting a high pressure region in the 
center follow behind. 

Computations are carried out with CFL — 0.8 and the 
Van Leer limiter, eq. 1651 . Spurious oscillations in vicinity of 
strong shocks are reduced by switching to the more diffusive 
minmod limiter, see ilAll No contact waves are present in 
the problem and, not surprisingly, the quality of our solution 
is essentially the same obtained by previous authors: the fast 
shocks are resolved in 2 -j- 3 cells, whereas the slow shocks 
are smeared out over 5-^6 zones. Very similar patterns are 
observed in the work of BA and dZBL. 

It is well known that Godunov-type schemes suffer from 
a common pathology, often found in these type of problems. 
In the classical cas e, this has been recognized for the first 
time bv lNohl jl987f) . The wall heating problem, in fact, con- 
sists in an undesired entropy buildup in a few zones around 
the point of symmetry. Our scheme is obviously no excep- 
tion as it can be inferred by inspecting the density profile in 

Fig.m 

We repeated the test with the HLL scheme and found 
that this pathology is worse when the HLLC scheme is used. 
The relative numerical undershoot in density, in fact, were 
found to be ~ 5% for the HLL and ~ 12%' for the HLLC 
scheme. Since similar errors were also reported by BA, and 
the same conclusions have been drawn in paper I, we raise 
the question as to whether the degree of this pathology grows 
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Figure 8. Rclativistic shock reflection problem at t = 0.4 on 1600 computational cells. The initial Lorcntz factor is 7 ta 22.4. Integration 
has been carried with the Van Leer limiter (except near strong shocks where the minmod limiter was used) and a Courant number of 
0.8. Notice the wall heating problem, evident in the density profile. 
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with the complexity of the Riemann solver. Future, more 
specific works should address this problem. 



4.3 Two-dimensional test problems 

Multi-dimensional numerical computations of magnetized 
flows are notoriously more challenging, due to the neces- 
sity to preserve the divergence-free constraint JjJ. In what 
follows, we consider three test problems: a cylindrical blast 
wave test, the interaction of a strong magnetosonic shock 
with a cloud and the propagation of an axisymmetric jet in 
cylindrical coordinates. 



4-3.1 Cylindrical Blast Wave 



Cylindrical explosions in cartesian coordinates are partic- 
ular useful in checking the robustness of the code and 
the algorithm response to different kinds of degeneracies. 
Here we follow the same setup adopted by KO, where the 
square [—6,6] x [—6,6] is filled with a uniform (p = 10 -4 , 
p g = 3 ■ 10~ 5 ), initially static (v = 0) medium, threaded by 
a constant magnetic field B = (B x ,0). The circular region 
\J x 2 + y 2 < 0.08 is initialized with constant higher density 
and pressure values, p = 0.01 and p g — 1 decreasing linearly 
for 0.08 r ^ 1. We adopt the ideal equation of state l|15[l 
with specific heat ratio F = 4/3. We consider two setups, 
corresponding to a relatively weak magnetic field B x = 0.1 
and a strong field B x — 1. Figures l9l and Hoi show the mag- 
netic field distribution, thermal pressure and Lorentz fac- 
tor for the two configurations at t = 4. Computations are 
carried using the van Leer limiter, eq. (1651 . together with 
the multidimensional limiting procedure described in i|A2l 
on 200 x 200 uniform grid zones. The Courant number is 
0.4. 

The expanding region is delimited by a fast forward 
shock propagating (nearly) radially at almost the speed of 
light. In the weak field case, a reverse shock delimits the 
inner region where expansion takes place radially. Magnetic 
field lines are squeezed in the y direction building up a shell 
of higher magnetic pressure. In the x direction the motion 
of the gas is not hindered by the presence of the field and it 
achieves a higher Lorentz factor (7 max = 4.39). In the strong 
field case, the expansion is magnetically confined along the 
x direction and the outer fast shock has reduced amplitude. 
The maximum Lorentz factor is 7 max = 4.02. 

We point out that numerical integrations for this test 
were possible only by locally redefining the total energy at 
the end of the time step: 



E -» E + 



'fa 



(76) 



where B c is the cell-centered magnetic field obtained after 
the Godunov step, whereas B f a is the new magnetic field 
obtained by averaging the face centered values given by (Eg. 
Notice that equation 1761 1 only redefines the energy contri- 
bution of the magnetic field that is not directly coupled to 
the velocity, see eq. 11211 and thus may be regarded as a 
first-order correction. In this respect, the energy correction 
we propose is the s ame u s ually adopted in CT schemes, see 
iBalsara fc Spicer I (1199911 . iTothl Jl997t) . Although this op- 
tional step results in a slight loss of energy conservation at 



Figure 11. Density gray scale map of the interaction between a 
strong shock and a cloud at t = 1. The upper and lower halves 
show the solutions computed with HLLC and HLL solvers, respec- 
tively, on 400 X 200 zones, with CFL = 0.4 and the MC limiter. 
Shock-flattening has been used to prevent spurious oscillations in 
proximity of the slow moving shock. 



the discretization level, it was nevertheless found to become 
particularly useful in problems where the magnetic pressure 
dominates over the thermal pressure by more than two order 
of magnitudes. 



4-3.2 Relativistic Shock-Cloud Interaction 

The interaction of a strong relativistic fast shock with 
a cloud is considered on the unit square [0, 1] x [0, 1] 
in 2-D cartesian coordinates (x,y). This problem has 
been extensively used for t e sting classic al MHD codes 
see (iDai fc Woodward 1 11994 ITothl Il99l and references 
therein). Here we consider a relativistic extension adopting 
a somewhat different initial condition, with magnetic field 
orthogonal to the slab plane. The shock wave travels in the 
positive x-direction and is initially located at x = 0.6. Up- 
stream, for x > 0.6, the flow is highly supersonic with pre- 
shock values given by (p, 7x,p 9 , -B z )pre = (1, 10, 10~ 3 ,0.5), 
where ^ x — (1 — t^) - ^. In this reference frame, shocked 
material is at rest with values given by 




/ 42.5942 \ 
127.9483 



post 



V 



-2.12971 



(77) 



Notice that the magnetic field carries a rotational disconti- 
nuity and the compression factor of density across the shock 
in not limited to 7 (we use T = 4/3) as in the classical case, 
but achieves a much higher value (~ 43). This feature is 
unique to relativistic flows. 

A circular density clump with p — 10 and radius 
r = 0.15 is placed ahead of the shock front, centered at 
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Figure 9. Gray scale levels of the x component of magnetic field (top left), y component of magnetic field (top right), gas pressure 
logarithm (bottom left) and Lorentz factor (bottom right) for the cylindrical blast wave with relatively weak magnetic field at t = 4. 
Magnetic field lines are plotted on top of the Lorentz factor distribution. Following KO, we use 32 equally spaced contour levels between 
0.008 and 0.35 (for B x ), -0.18 and 0.18 (for By), -4.5 and -1.5 (for Logp 9 ), 1 and 4.57 (for 7). 



(x,y) — (0.8,0.5). Transverse velocities v y and v z and the 
x and y components of magnetic field are set to zero every- 
where. We use 400 x 200 computational zones, by assuming 
reflecting boundary at y — 0.5 and free flow across the re- 
maining boundaries. The MC limiter, eq. (1631 . is employed 
everywhere except in proximity of strong shocks where we 
revert to the minmod limiter, see ^1 All The Courant number 
is 0.4. 



Shortly after the impact, the cloud undergoes strong 
compression with the density rising by a factor of more than 
20. The collision generates a bow fast shock propagating in 
the shocked material and a reverse shock is transmitted into 
the cloud. After the transmitted shock reaches the back of 
the cloud, the two bent parts of the original incident shock 
join back together and complicated wave pattern emerges. 
By t — 1 the cloud is completely wrapped by the incident 
shock, and the cloud expands in the form of a mushroom- 
shaped shell, see upper half of Fig. 1111 The solution com- 
puted with the HLL solver (lower half in Fig. II IP show sim- 
ilar structures, although the amount of numerical viscosity 
is considerably higher. 



Notice that, because of the assumed slab symmetry, the 
condition v ■ B = is preserved in time and the solution to 
the Riemann problem at each interface consists of a three 
wave pattern: two fast waves separated by a tangential dis- 
continuity. In this regard, our HLLC solver provides a better 
approximation of the full wave structure. 



4-3.3 Relativistic Jet 

As a final example, we consider the propagation of an ax- 
isymmetric jet in cylindrical coordinates (r, z). The con- 
figuration adopted he re corresponds to model C2-pol-l in 
iLeismann et alJ (I2005T) . 

The domain [0, 12] x [0, 50] (in units of jet beam) is 
initially filled with a static uniform distributions of density, 
gas pressure and magnetic field, given respectively by 



Pa 



Pa 



r(r- i)M 2 



•r»r 



B z 



(78) 



The numerical value of p a follows from the definitions of the 
beam Mach number M — Vb/c 3 = 6, jet to ambient den- 
sity ratio 77 = 10 -2 and beam axial velocity Vb = 0.99. 



An HLLC Solver for Relativistic Flows 15 




The ideal equation of state 1151 is used with V = 5/3. 
The jet nozzle is located at the lower boundary r ^ 1, 
z — 0, where boundary conditions are held constant in time, 
(p,v r ,v z ,B r ,B z ,p g ) = (rj,0,v b ,0, B z ,p a ). For r > 1 we pre- 
scribe boundary values with antisymmetric profiles for ax- 
ial velocity and radial magnetic field. Symmetric profiles are 
imposed on the remaining quantities. This configuration cor- 
responds to a twin counter jet propagating in the opposite 
direction. Outflow boundaries are imposed on all other sides, 
except at r = where reflecting boundary conditions are 
used. We employ a uniform resolution of 20 zones per beam 
radius and carry integration until t — 1 26 with CFL — 0.4. 

The results are shown in Fig. 1121 where we display 
density logarithm (upper panel), magnetic pressure (mid- 
dle panel) and Lorentz factor distributions (lower panel). In 
each panel, the upper and lower halves show the solutions 
obtained with the HLLC and HLL solvers, respectively. As 
we already pointed out in the non magnetic case (Paper 
I), the HLLC integration features considerably less amount 
numerical diffusion as evident from the richness in small 
scale structures, notably in the density distribution. In fact, 
density is the physical quantity more sensitive to the in- 
troduction of the tangential wave i n the Riemann solver . 
Comparing our results with those of llLeismann et aljl2005l 



see their Fig. 5) we can observe that our solution has a sim- 
ilar (or even larger) richness in fine structure details at half 
the resolution (20 ppb in our case, 40 ppb in their case). 



5 CONCLUSIONS 

An HLLC approximate Riemann solver has been developed 
for the relativistic magnetohydrodynamic equations. The 
new approach improves over the single state HLL solver in 
the ability to capture exactly isolated tangential and con- 
tact discontinuities. Several test problems in one and two 
dimensions demonstrate better resolution properties and a 
reduced amount of the numerical diffusion inherent to the 
averaging process of the single state HLL scheme. The solver 
is well-behaved for strictly two-dimensional flows, although 
applications to genuinely three-dimensional problems may 
suffer from a pathological singularity when the component 
of magnetic field normal to a zone interface approaches zero. 
This feature does not persist in the classical limit. 

Multidimensional integration has been formulated in 
a versatile and efficient way within the framework of the 
corner transport upwind (CTU) method. The algorithm 
is stable up to Courant numbers of 1 and preserves the 
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Figure 12. Gray scale images of density (top panel), magnetic 
pressure (middle panel) and Lorentz factor (bottom panel) for 
the axisymmetric jet. The upper (lower) half in each panel refers 
to the integration carried with the HLLC (HLL) solver. Both 
integrations were carried till t = 126 with CFL = 0.8 and the 
Van Leer limiter. An ideal equation of state is used with Y = 5/3. 
Magnetic field lines are plotted on top of the Lorentz factor gray- 
scale images. 

divergence-free condition via constrained transport evolu- 
tion of the magnetic field. The additional computational cost 
and the numerical implementation in an existing relativistic 
MHD code are minimal. 
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APPENDIX A: 

Al Shock Flattening 

For strong shocks, we found that the one-dimensional pre- 
scriptions 163H or 165H can still produce spurious numerical 
oscillations eventually leading to the occurrence of negative 
pressures. A weak form of flattening is introduced by re- 
placing eq. 16311 or 1651 with the minmod limiter whenever 
a strong shock is detected. In order for the latter condition 
to hold, we require that both V • v < and Xmtn = 0, where 
V ■ v is computed by central differences whereas 

Xmin = mm(xf + ij,X?,j>Xi : -i,j)Xi,j + i>Xij»Xi,j_i) ■ ( A1 ) 
The switches \ x an d X v are designed as follows 

Pi+1,3 -Pi-l.j < 

min(p<+i,j,pi_i,j) ^ ' (A2) 
otherwise , 

Pij+i -Pi,j-i < 
min(p ij+ i,pij_i) " ( A3 ) 

otherwise , 

where we set e = 5 in all computations presented in this 
paper. 




A2 Multidimensional Limiting 

Occasionally, we found that strong shocks propagating 
obliquely to the grid in highly magnetized media may ben- 
efit from an additional form of limiting, based on genuinely 
multidimensional constraints. When needed, we enforce the 
maximum and minimum interpolated values in each cell 
(i, j) to lie within the bounds provided by the four neighbor- 
ing zones (i + l,j), (i - + 1), - 1). Specifically, 
denote with <J max and <J mm the maximum and minimum val- 
ues of q £ V in these cells. Once the limited slopes 8 x q and 
5 y q have been computed according to 163H or 1651 . we apply 
the correction 

S x q — * r5 x q , S v q — » rS y q , (A4) 

where the mu lti-dimensional limiter r is constructed as in 
iBalsaral 12004) : 

^ = ^(l^^(^^-Sf-)) ' (A5) 

with 5 max = WBx(\5 x q\, \5 v q\), 5 min = mm(\S x q\,\8 y q\). We 
set ip = 2 for density and magnetic field, ip = 3/4 for velocity 
and i/j = 1 for thermal pressure. 



